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ABSTRACT 


This  paper  describes  a  new  statistical  approach  to  image  segmentation. 
Making  nse  of  Gibbs  distribution  models  of  Markov  random  fields  a  dynamic 
programming  based  segmentation  algorithm  is  developed.  A  number  of 
examples  are  presented  which  give  an  indication  of  the  potential  of  this 
approach. 


I. 


INTRODUCTION 


3  This  report  presents  a  new  statistical  approach  to  the  image 
segmentation  problem.  By  modelling  image  data  as  a  Markov  random  field 
characterised  by  a  Gibbs  distribution,  a  dynamic  programming  algorithm  is 
developed.  The  primary  contribution  of  the  paper  is  this  new  near  optimal 
method  for  processing  scenes  described  by  the  non-causal  Gibbs  model.  ^ 

Image  segmentation,  the  process  of  grouping  image  data  into  regions 
with  similar  features  is  a  component  process  in  image  understanding  systems 
and  also  serves  as  a  tool  for  image  enhancement.  As  such  it  has  received 
considerable  attention  in  the  literature.  Many  techniques  work  well  on 
noise  free  images  with  slow  spatial  variation  in  intensity.  However,  when 
the  data  is  noisy  or  textured,  these  algorithms  become  less  reliable.  In 
this  case  it  can  be  advantageous  to  statistically  model  the  noise  and  any 
texture  which  is  random  in  nature.  Furthermore,  one  must  also  take 
advantage  of  two-dimensional  spatial  ergoticity  to  average  the  effects  of 
noise.  If  a  region  is  spatially  ergodic  then  a  pixel  and  its  neighbors 
will  have  similar  statistical  properties.  In  its  simpler  forms,  the  Gibbs 
model  can  be  used  to  exploit  this  type  of  spatial  continuity,  and  this  is 
its  primary  role  in  the  segmentation  algorithm. 

Use  of  the  Gibbs  distribution  dates  back  to  the  work  of  Ising  [1]  in 
1925  who  modelled  molecular  interaction  in  ferromagnetic  materials,  and  it 
has  received  considerable  attention  in  both  the  statistical  mechanics  and 
statistics  literature  [2].  However,  only  recently  have  attempts  been  made 
to  apply  it  to  problems  in  image  processing.  In  [4],  the  autobinomial  form 
of  the  Gibbs  distribution  was  used  to  model  texture.  The  algorithm  in  [S] 
segments  textured  images  hierarchically  operating  on  successively  smaller 
blocks  and  uses  Gibbs  distributions  to  model  texture.  To  our  knowledge, 
the  work  in  [6]  represents  the  first  application  of  the  Gibbs  model  to 
Image  Segmentation.  The  algorithm  in  [6]  is  highly  parallel  in  nature  with 
the  flavor  of  a  'relaxation'  algorithm  and  requires  a  number  of  iterative 
passes  on  the  image  data.  The  algorithm  we  propose  processes  the  data  in  a 
raster  scan  fashion,  and  only  requires  a  single  scan  of  the  data.  It  will 
be  iaqtortant  to  further  study  the  trade-offs  in  the  various  algorithms. 

The  report  is  organized  as  follows.  Section defines  the 
segmentation  problem  in  a  statistical  framework,  introduces  some  notation, 
and  presents  some  background  on  Markov  random  fields.  Section-fTT  then 
presents  the  dynamic  programming  algorithm  in  detail  for  the  case  of 
segmenting  images  consisting  of  uniform  intensity  regions  in  high  levels  of 
additive  white  Gaussian  noise.  Sect ion-fT* presents  results  of  applying  the 
algorithms  to  some  experimentally  generated  images  consistent  with  this 
model  as  well  as  some  synthetic  aperture  radar  images  which  are  clearly 
inconsistent  with  the  assumed  model.  These  results  clearly  demonstrate  the 


applicability  of  the  technique  to  realistic  data  as  well  as  the  robustness 
of  the  algorithm  with  'spect  to  modelling  assumptions.  In  Section^  some 
comments  and  concluding  remarks  are  given,  and  extensions  to  this  work 
which  are  in  progress  are  1  iefly  outlined. 


s. 


II.  PROBLEM  FORMJLA'pION  AND  MATHEMATICAL  BACKGROUND 


Preliminary  Definitions 


Let  *  class  of  scenes  be  characterized  by  a  discrete  finite  random 
field  X  =  [X.jJ  of  size  (N^XN^),  and  let  a  realization  of  this  field  or  a 

specific  scene  be  represented  by  the  matrix  x  =  [x.  J.  It  will  be  assumed 

that  each  pixel  (i»  j)  can  belong  to  one  of  M  distinct  region  types  and 
that  X^  a  a  if  pixel  (i,  j)  is  a  member  of  region  m.  me[l,  2,  ....  M] . 

Associated  with  a  specific  scene  is  a  set  of  K,  (N^XN^)  observation 


matrices  y  =  {yk  }  *  ,  yk 


[y*  ].  For  simplicity  of  exposition,  we  will 


assume  K  ■  1  and  simply  define  y  to  be  y  =  [ y4 j 1  *  However,  it  should  be 

pointed  out  that  the  algorithms  presented  below  extend  trivially  to  the 
case  of  multiple  observations  snch  as  with  Landsat  data.  Since  regions  can 
be  textured  or  contain  observation  noise  the  range  space  of  y^  is  larger 

than  that  of  X^.  Thus  y  will  be  assumed  to  be  a  realization  of  a  real 

valued  random  field  T  =  [Y^l.  The  general  model  which  can  we  will  employ 

ia 


Yii  *  VV  +  wu 


(1) 


The  field  V  is  a  random  noise  field,  and  the  mapping  F..  can  be  used 
•  J  ^  J 

to  characterize  texture  models.  In  the  case  which  will  be  described  in 
most  detail,  a  region  will  be  characterized  by  constant  intensity  so  that 


r 

f: 

Jr. 


VV  ■ -  * 


<2> 


Furthermore,  N^  will  be  assumed  to  be  a  white  Gaussian  field  with  zero 


mean  and  variance  a  ,  i.e. 

W. 


'ij  -^/(0 


.  .2> 


(3) 


The  segmentation  problem  can  nov  be  simply  stated  as  follows.  Given 
the  observation  matrix  y  =  fy^]  find  an  estimate  x  =  [x^]  of  the  scene 
realisation  x  =  [x.^J.  The  aljorithm  presented  below  attempts  to  maximize 

the  posterior  probability  or  likelihood  of  x  given  y.  In  particular  if 
P(.)  is  an  appropriate  probability  measure,  then  one  would  like  to  find  the 

estimate  x  which  maximizes  P(X  =  x  I  Y  =  y) ,  Using  Bayes  rule 

P(X=I  |  Y=y)  =  P(Y=y  I  X=x)  P(X=x)  (4) 

P(Y=y) 

Since  P(Y-y)  is  independent  of  the  estimate  x,  we  can  equivalently  maximize 

P(X~;.  Y=y)  =  P(Y=y  I  X  =  i)  P  (X  =  I)  (5) 

or 

In  P  (X-i.  Y=y)  =  In  P(Y=y  |  X=J) 


+  In  P(X=£)  (6) 

The  dynamic  programming  algorithm  presented  in  the  next  section  is  an 

approximation  to  one  which  guarantees  finding  the  x  which  maximizes  (6). 
It  should  be  pointed  out  that  the  difficulty  in  maximizing  (6)  is  that  it 
is  a  joint  log-1 ike 1 i hood  for  all  the  image  data.  It  does  not  simply 
describe  the  likelihood  of  a  single  pixel.  In  particular  the  maximizing 

N1N2 

x,  is  one  of  N  possibilities. 


and  the  Gibbs  Distrlbut ion 


Obviously,  any  processing  algorithm  for  maximization  of  (6)  will 
depend  critically  on  the  form  of  P(X=x)  and  P(Y=y  I  X=x).  In  this 
subsection,  these  measures  are  defined  in  the  case  where  X  is  a  Markov 
random  field  characterized  by  a  Gibbs  distribution,  and  Y  is  given  by  (1)  - 
(3). 
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To  begin  it  «i;1  be  helpful  to  introduce  some  additional  notation. 

Let  L  de  defined  as  .  -  N^XN,,  lattice  characterizing  all  pixel  locations  in 

the  scene,  i.e. 

L  =  { (i, j) :  1  <  i  <  Nj.  1<  j  <  N2)  (7) 

Next  let  define  a  set  of  neighbor  sites  for  the  pixel  (i,  j)  but 

excluding  (i,  j)  itself.  Tiro  siaple  neighborhoods  are  depicted  in  figure 
1.  These  are: 

=  ( ( 1 , a) :  0  <  <i  -  l)2  +  (j  -  a)2  <  1}  (8) 

Hj2  -  {<1,  a):  0  <  (i  -  l)2  +  (j  -  a)2  <  2}  (9) 

Finally,  define  to  be  the  collection  of  all  neighborhoods  in  L,  or  a 

1  2 

neighborhood  systea.  Then  ^  and  q^  characterize  the  collection  of  all 
1  2 

neighborhoods  q^  and  q^  ,  i.e. 

nL  =  {,,ij  :  (1,j)  8  L}  (10) 

q2  *  ^ij  :  8  (ID 

Given  this  notation,  a  Markov  randoa  field  can  then  be  defined  as 


A  randoa  field  on  a  lattice  L  is  a  aarkov  randoa  field  with  respect  to 
a  neighborhood  systea  q^  if  and  only  if 


P(Xij  =  xij  1  Xla  =  Xla*  (1,  a)  8  L’  (1,B) 
■  P(Xij*xiJ  1  xu‘  *1.- 


If  in  addition  to  (12)  P(X  =  x)  >  0  for  all  realizations  x,  then  X  can  be 
characterized  by  a  Gibbs  distribution  defined  on  the  neighborhood  systea 
HL(2l ,  [3].  In  order  to  define  the  Gibbs  distribution  it  is  necessary  to 

introduce  the  notion  of  the  cliques  of  a  neighborhood  system.  Simply,  a 


8 


As  can  be  seen  fros  the  definition  Z  is  simply  a  normalising  constant  so 
that  the  sum  of  the  probabilities  of  all  realizations,  x,  add  to  one.  Thus 
the  key  functions  in  determining  the  properties  of  the  distribution  are  the 
potential  functions  Vc(x).  The  only  limitation  on  Vc<x)  is  that  it  only 

depend  on  the  values  of  the  pixels  in  clique  c.  Ve  will  consider 

homogeneous  fields  where  the  form  of  V  (x)  is  fixed  by  the  structure  of  the 

c 

clique  c  and  not  its  location  in  the  lattice.  For  the  segmentation 
algorithm  discussed  in  the  next  section,  the  potential  functions  were 
chosen  to  exploit  spatial  continuity,  and  for  simplicity  of  calculation. 

Ve  have  modelled  the  region  clustering  as  being  characterized  by  a  Markov 

2 

random  field  over  the  neighborhood  system  q^.  The  most  general  form  we  use 

for  the  potential  functions  V  (x)  are  given  in  Table  1. 

c 


Cl iqne 

No.  Structure  Potential 

1  •  I(i.j)]  V  (x)  =  a  if  x..=  m 

c  n  ij 

2  **  I(i» j) »(i,  j+1)]  V  (x>*  p.  if  x.  =  x. 

®  ij  >  J  +  l 

-pj  otherwise 

3  •  ( (i> j)  .  ( i-1 , j) 1  V  (x)  =  p,  if  x.  =  x.  ,  . 

c  2  ij  i  +  1,  j 

-?2  otherwise 

4  .•  [<i.j).<i-l.j+l)]  V  (x)  =  p.  if  x..  =  x.  ,  _ 

c  3  ij  i-1, j+1 

-pj  otherwise 

5  %  Ki.j).<i+1,  j+1)]  V  (x)  =  p.  if  x.  =  x.  A  , 

«  4  ij  l  +  1,  j  + 

~P^  otherwise 

6  **  Ki.j),  (i+1*  j).  Vc(x>  m  if  all  x^  in  c  are 

(i»  j+1)]  equal 

“7^  otherwise 

7  Ki.j),( i-1 , j) ,  (  i. j+1) ]  V  (x)  =  r,  if  all  x.,  in  c  are 

*  •  ij 

equal 

~7j  otherwise 

*  **  Ki.j)  ,(i,  j+1)  , (i+1, j+1)]  V  (x)  =  7,  if  all  x  in  c  are 

c  a  IJ 

equal 

~7j  otherwise 

9  ( (i, j) , (i, j+1) , ( i-1, j+1) ]  V  (x)  =  7.  if  all  x4J  in  c  are 

C  4  ij 

equal 

~7^  otherwise 


10  ::  Ki.j) , (i-i, j) . c i, j+D, 

( i, j+1)] 


V  (x)  =  a.  if  all  x.  in  c  are 
c  1  ij 

equal 


-Oj  otherwise 

T  LEI 

Potential  Funct-c^s  for  Seguentation  Algorithm 


The  parameters  a^,  p.,  y  o,  need  to  be  estimated  for  certain  classes 

of  scenes.  Although  this  is  a  difficult  problem,  there  are  methods  in  the 
literature  for  estimating  these  parameters,  see  [3]  -  [5].  However,  we  are 
not  using  the  Gibbs  distribution  to  model  textures  or  detailed  shape.  We 
just  want  to  model  the  fact  that  regions  are  clusters  of  pixels,  i.e. 
spatial  continuity  within  regions.  For  the  examples  in  section  IV,  we  use 
only  cliques  1-5,  an*,  in  many  cases  just  cliques  2  and  3,  e.g.  y .  =o .  =  0. 

The  nonzero  parameters  were  chosen  by  trial  and  error.  Although  for  a 
given  signal  to  noise  rati  the  algorithm  was  relatively  insensitive  to  the 
choice  of  these  values,  as  the  signal  to  noise  ratio  changed  we  found  it 
necessary  to  modify  the  values  We  are  presently  studying  this  phenomena 
more  carefully  as  well  as  determining  new  schemes  for  estimating  these 
parameters  for  the  context  in  which  they  are  being  used. 


Finally,  since  a  particular  realization  of  the  Gibbs  field,  x,  assigns 
each  pixel  to  one  of  the  M  region  types,  using  the  Gaussian  noise  model  an 
appropriate  form  for  P(Y=y  |  X=x )  is 


P(Y=y  IX-x)  =TT 
m=l 


(  i, j) eS 


(2no*)x/* 


exp(z— r  (y.  -r 
2c*  -7ij  i 


)2) 


Sm=  1  Xij  "■  } 


(16) 

(17) 
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III.  NEAR  OPTIMA.  tfiP  SEGMENTATION 


In  this  section,  an  ptimal  algorithm  is  posed  for  processing  images 
consisting  of  a  few  rows.  The  complete  near-optimal  algorithm  is  then 
obtained  by  applying  this  opt;mal  processor  on  overlapping  strips  of  the 
larger  NjXN^  This  algorithm  will  be  near  optimal  when  correlation 

between  the  random  variables  in  X^  drops  rapidly  as  their  vertical 

distance  increases.  This  assumption  appears  reasonable  in  the  sense  that 
it  can  often  be  shown  for  one  dimensional  Markov  chains.  However, 
calculation  of  correlations  in  a  Markov  random  field  is  extremely  difficult 
even  for  the  simplest  case  (Ising  Model  [2]).  and  is  an  unresolved  problem 
in  the  statistics  literature. 

First  consider  the  problem  of  segmenting  a  DXN2  image,  D  <<  .  Let 

£(.)  denote  a  log-1 ike 1 i hood  function.  Using  the  Gibbs  likelihood  and  the 
conditional  data  likelihoods  derived  in  the  last  subsection  of  section  II, 
we  can  write  the  following  expression  for  the  joint  log-1 ikel ihood  of  the 
image  data  y  and  a  realization  x  of  the  Gibbs  field  X: 


UY- 


y.  X=x)  =  £(Y=y  I  X=x)  +  £(X=x) 


S,(X=x)  =  -InZ  +  ^ 


V  (x) 


(18) 


(19) 


ceCL(ti* ) 


DN„ 


M 


t<r-,  I  1-r)  -  -y3  i.(2,„2>-  }  }  £  ‘y.j-V2 

m=l  (i,j)eS 


(20) 


Observe  that  we  can  calculate  (18)  recursively  as  follows 

iLj  =  «.( Y=y,  X=x ) 

1 


-DN 


l  = 


—  ln(2no^)  -In  Z 


M 


•  *1-1  +  2  V*>  }  2^ 


v_l  2 

C  '  =  { c sC v  ):  c  contains  only  pixels  in  column  k  or  only  pixels 

in  columns  k  -  1  and  k  ) 


S  =  {  (i.j):  x. .  =  m.  j  =  k)  1  <  i  <  f 
m  »J 

This  recursion  in  conjunction  with  the  principle  of  optimality  allows 

formulation  of  a  forward  dynamic  programming  algorithm  [7]  for  finding  x. 
The  state  space  associated  with  the  dynamic  programming  algorithm  has 

dimension  since  there  are  possible  segmentations  of  each  column  of 
the  DXN^  scene.  This  implies  that  the  algorithm  would  have  Nj  iterations 

with  on  the  order  of  M20  calculations  during  each  iteration.  Thus 

this  algorithm  may  only  be  computationally  tractable  for  small  values  of  M 
and  D,  e.g.  2  <  M,  D  <.  4.  Although  using  the  technique  described  below, 
full  size  images  can  be  processed  with  reasonable  computation  speeds,  we  do 
only  recommend  the  algorithm  be  implemented  for  segmenting  images  into  at 
most  MM  region  types.  Ve  are  preaently  working  on  a  method  which  allows 
this  algorithm  to  be  applied  iteratively  to  segment  images  for  which  M  >  4. 
More  will  be  said  regarding  this  in  Section  V.  Since  this  is  a  standard 
dynamic  programing  application,  details  will  not  be  given.  However,  two 
important  remarks  are  in  order. 

Remark  1 


Observe  that  the  value  of  o  is  independent  of  any  segmentation  x,  and 

hence  the  algorithm  can  be  initialized  by  setting  o  =  0.  In  particular, 

there  is  no  need  to  undertake  the  difficult  task  of  calculating  the 
partition  function  Z. 


Remark  2 


In  order  to  calculate  potentials  for  all  cliques  in  column  one  and  row 
one  of  the  (DX^)  image,  it  is  necessary  to  assume  a  segmentation  for  a 

ficticious  column  zero  and  row  zero.  This  corresponds  to  a  boundary 
condition  for  the  Markov  random  field  X.  An  appropriate  choice  for  these 
boundary  conditions  will  be  discussed  below. 


In  order  to  use  the  dynamic  programming  algorithm  described  above 
which  is  capable  of  optimally  processing  a  DXN^  strip  of  an  NjXN^  image 

N^>>D,  ve  will  assume  the  random  variables  X^  and  X^+p  ^  for  all  (i,j)  to 

have  negligible  correlation  or  covariance.  Thus,  the  segmentation  for  row 
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A 


i  should  be  negligib  impacted  by  the  segmentation  of  row  i+D.  In  view  of 
this,  the  complete  se^,  ntation  procedure  is  described  by  the  following 
algorithm. 


To  summarize,  the  dynamic  programming  algorithm  is  applied  to 
overlapping  image  strips  of  width  D,  but  only  the  segmentation  of  the  first 
row  of  that  strip  is  used.  For  example,  the  processing  of  rows  1  through  D 
yields  a  segmentation  for  row  1,  and  the  processing  of  rows  2  through  D+l 
yields  a  segmentation  of  row  2.  Under  the  correlation  assumption  above, 
this  algorithm  is  near  optimal  since  the  data  in  row  I+D  will  have  little 
impact  on  the  segmentation  of  row  I. 


As  pointed  out  in  Remark  2  above,  the  strip  dynamic  programming 
algorithm  requires  a  fixed  segmentation  or  boundary  condition  for  row  1-1 
and  column  0.  For  1=1,  (i.e.  row  0)  and  column  zero,  we  arbitrarily 

where  N^e  [1, 

2,  ...,  N]  and  is  the  assumed  background  intensity.  Although  this  is 

somewhat  arbitrary,  if  the  correlation  assumption  holds  both  vertically  and 
horizontally,  it  will  have  negligible  impact  on  the  segmentation  of  the 
scene  for  i  >  D  and  j  >  D.  For  I  >  1  we  use  the  fixed  segmentation  of  the 
previous  row  to  initialize  the  strip  processor. 


assume  all  pixels  to  be  backround  pixels,  i.e.  x^  _  x^  = 


u 


IV.  EXAMPLES 


In  this  section,  some  examples  are  presented  which  are  representative 
of  the  performance  of  the  algorithm  and  which  highlight  some  of  its 
properties.  To  begin,  conside:  figures  3  and  4.  They  show  the  results  of 
applying  the  algoirthm  to  (64  X  64)  test  images  consisting  of  an  object 
either  an  ellipse  or  a  diamond,  on  a  background.  Thus,  M=2 .  For  these 
examples,  the  images  were  corrupted  by  additive  white  zero  mean  Gaussian 
noise  fields  with  variances  such  that  the  signal  to  noise  ratio  S/N=2  where 
we  define 


a 


The  values  of  r^,  r^,  and  a  were  assumed  known  and  the  algorithm  was 

applied  with  D=4  and  all  Gibbs  parameters  set  to  zero  except  the  p..  For 

the  ellipse  of  Figure  3  P .  =  0.2  while  for  the  diamond  of  Figure  4 

p.  =  0.15  i=1.2»3,4.  The  top  row  of  each  figure  from  left  to  right  show 

the  original  and  noise  corrupted  images,  while  the  bottom  row  from  left  to 
right  show  the  Gibbs  algorithm  segmentation,  and  the  result  of  filtering 
the  segmented  image  through  a  (3  X  3)  median  filter.  As  can  be  seen  by 
comparison  of  the  two  figures,  the  algorithm  performs  well  at  this  signal 
to  noise  ratio  but  as  the  magnitude  of  the  Gibbs  parameters  decreases,  the 
algorithm  becomes  subject  to  more  single  pixel  errors.  This  is  consistent 
with  expectation  since  in  this  case  the  Gibbs  distribution  is  weighted  less 
relative  to  the  data  terms  in  the  likelihood  function.  Thus  there  is  less 
emphasis  on  spatial  continuity  in  the  likelihood. 

Figures  5-7  show  the  result  of  applying  the  algorithm  to  an  ellipse 
in  noise  such  that  S/N  =  1  using  D«2,3,4.  Although  performance  is  best  for 
DM  there  is  surprisingly  little  degredation  in  performance  as  D  decreases. 
In  all  three  figures  p. =0.175,  i=l,2,3,4  and  all  other  Gibbs  parameters 

were  zero. 

Figures  8  and  9  show  the  results  of  applying  the  algorithm  to  a 
diamond  in  noise  such  that  S/N=l,  In  both  figures  DM ,  however  by 
comparison  of  figures  8  and  9,  one  can  see  the  improvement  in  performance 
which  can  be  obtained  by  assuming  some  additional  prior  knowledge  of  shape. 
For  Figure  8  p^O.2  i-1,2,3,4  and  all  other  Gibbs  parameters  were  zero. 

For  Figure  9,  realizing  that  the  object  had  diagonal  edges,  the  diagonal 
cliques  were  emphasized  relative  to  the  horizontal  and  verticle  cliques, 
i.e.  p  • 1 5  »»d  Pj»p^“0.25  while  all  other  parameters  were  zero.  For 


the  letter  case,  !ht  shape  of  the  diamond  was  more  obvious  in  the 
segmentation. 


For  Figures  10  -  17,  he  algorithm  was  applied  to  two  test  images, 
each  containing  B=4  region  types.  For  Figures  10  -  IS  D=2  while  for 
Figures  16  and  17  D=3 .  Ve  usel  D=2  to  keep  computation  times  down,  however 
as  can  be  seen  from  Figures  12  -  IS,  while  performance  at  S/N=1.5  is 
reasonable,  performance  at  S/h^l  is  not.  For  M=4,  we  define  the  signal  to 
noise  ration  as 


S/N 


Increasing  to  D=3  did  considerably  improve  performance  at  S/N=l,  however 
CPU  time  for  this  case  on  a  VAX  11/780  was  60  minutes  as  compared  to  4 
minutes  for  16=4  and  D=2  and  7  minutes  for  N=2  and  D=4 .  This  is  our 
motivation  for  going  to  the  hierarchical  scheme  for  handling  multi-region 
images  (N  >  3)  which  will  be  briefly  outlined  in  Section  IV  below.  We  also 
point  out  that  we  are  rewriting  our  code  to  make  use  of  look-up  tables  and 
anticipate  a  factor  of  10  or  better  improvement  in  computation  speed. 


Finally,  for  Figures  19-21,  the  algorithm  was  applied  to  (64  X  64) 
sections  of  the  synthetic  aperture  radar  image  shown  in  Figure  18.  Figure 
19  shows  the  result  of  applying  algorithm  to  a  small  bay,  while  in  Figure 
20,  there  is  a  river  with  some  bridges,  and  in  Figure  21  there  is  a  boat. 

In  all  cases,  we  assumed  M=2 ,  and  obtained  r^,  r^  and  a  by  applying  the 

2 

method  of  moments  under  the  assumption  that  the  (64)  pixels  in  the  scene 
where  samples  of  a  random  variable  characterized  by  a  distribution  which 
was  a  mixture  of  two  Gaussian  distributions.  Each  picture  shows  the 
original  scene  and  three  segmentations  corresponding  to  different  Gibbs 
parameters.  For  Figure  19,  proceeding  clockwise  aronnd  the  picture  we  had 
<^=-.05  (white),  a^.OS  (black)  and  all  pJ=0.125,  ^=<*2=0  and  all  p. =0.125, 

*nd  all  0j=O.15.  For  Figure  20,  we  had.  <»2=.05  (white),  U2=-.05 
(black),  and  all  0.=O.1,  Oj=.05,  o^-.OS,  »nd  all  0^=0. 125,  ai=°2=®  and  all 
02=0.  For  Figure  21  we  had,  Oj=O2=0  and  all  02=0.2,  Oj=a2=0  and  all 
02=0.35,  Uj^O.l,  O2=~0.1  and  all  02=0.35. 


16 


V.  CONCLUDING  R-MA^KS 


This  report  has  prese  ted  a  new  approach  to  segmentation  of  noisy 
images.  It  uses  the  Gibbs  distribution  to  model  spatial  continuity  or 
clustering  properties  of  reg  iors.  The  algorithm  is  recursive  in  nature, 
requires  a  single  pass  on  the  data,  and  works  well  at  low  signal  to  noise 
ratios. 

Presently,  two  extensions  to  this  algorithm  are  being  developed.  The 
first  allows  computationally  efficient  segmentation  of  images  with  more 
than  two  region  types.  If  there  are  M  region  types,  the  two  region 
version  of  the  algorithm  is  applied  M-l  times.  Thus  computation  time  grows 
linearly  with  the  number  of  regions.  To  see  how  this  can  be  done  consider 
the  case  where  M=4  and  for  convenience,  the  regions  have  been  labelled  so 
that  r^<  Tj<  rj  <  r^.  The  two  region  algorithm  can  first  be  applied  using 

region  means  of  end  r ^ .  Since  r^  <  ^  and  r^  >  r^  the  result  will  be  a 

segmentation  grouping  regions  1  and  2  together,  and  regions  3  and  4 
together.  Next  the  two  region  algorithm  is  applied  to  the  image  using 
region  means  of  r^  and  rj.  However,  only  pixels  classified  as  being  in  the 

region  with  mean  r 2  during  the  first  pass  are  classified  in  the  second 

pass.  Those  classified  as  being  in  the  region  with  mean  r^  during  the 

first  pass  are  ignored.  This  yields  a  three  class  segmentation  into 
regions  1.  2  and  the  combined  regions  3  and  4.  To  separate  regions  3  and  4 
a  third  pass  is  made  using  region  means  r^  and  r^  but  ignoring  those  pixels 

assigned  to  regions  1  and  2  during  the  first  pass.  Ve  feel  that  this 
hierarchical  algorithm  will  have  little  effect  on  the  near  optimality  of 
the  overall  approach,  however  we  have  no  examples  to  show  at  this  time. 

The  second  extension  is  to  textured  images.  In  this  case,  instead  of 

aiodelling  regions  as  being  of  constant  intensity  but  imbedded  in 

uncorrelated  noise,  we  asstae  textured  regions  are  realizations  of  a  second 

Gibbs  distribution.  In  this  case,  we  are  developing  algorithms  similar  to 

the  one  given  above,  but  with  the  Gaussian  (quadratic)  data  term  in  the 

log-1 ike 1 i hood  replaced  by  terms  associated  with  the  second  Gibbs 

distribution  characterizing  texture.  In  this  context,  we  are  also 

developing  new  methods  for  Gibbs  parameter  estimation  in  texture  images. 

The  main  difference  between  our  work  on  texture  modelling  and  that  in  [4], 

[S]  is  that  we  are  experimenting  with  different  types  of  potential 

functions  V  (x) . 
c 

In  conclusion,  we  feel  that  the  applications  of  Gibbs  models  to 
problems  in  image  processing  and  image  understanding  are  just  beginning  to 
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Figure  4  -  Segmentation  of  a  diamond  with  S/N=2 


Figure  6  -  Segmentation  of  an  Figure  7  -  Segmentation  of  an  ellipse 

ellipse  with  S/N=l  and  0=3  with  S/N=l  and  D=4 


Figure  10  -  Segmentation  of  a  four  region  diamond  test 
Image  with  S/N=2  and  D=2 


Figure  12  -  Segmentation  of  a  four  region  diamond  test 
image  with  S/N=1.5  and  D=2 
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